function [means total]= default_setting_mult_wz(rho, alpha, delta,addz, z, n, j);
    % calculate the mean and total of simulated donations
    % inputs: 
    % rho: simulated from distribution of (nonzero) donation in AD treatment
    % lambda1, lambda2 and lambda3
    % addz: zeros added to represent ungenerous types.
    % simulate donations when default is set rho+j
    % alpha: simulated alpha costs
    % delta: simulated delta costs
    % z indicator for agents that faces deviation costs. 
    % n: sample size
    % j: set default = rho*j

def=rho.*j;

Vx= (rho.^2)/2;
Ux= Vx- delta;
U0= -alpha;
Udef= rho.* def- (def.^2)./2;

sim_donation=zeros(n,1);
% if no cost --> pick rho instead of default
sim_donation(~z)=rho(~z);
% if endure switching cost
% deviate if V(rho)-delta>V(d)
switchtorho=zeros(n,1);
switchtorho(z)= (Ux(z)>=U0(z) & Ux(z)>Udef(z));
% stay at default with cost
stay_default=zeros(n,1);
stay_default(z)= (Udef(z)>=Ux(z) & Udef(z)>= U0(z));

optout=zeros(n,1);
optout(z)=(U0(z)>Ux(z) & U0(z)>Udef(z));

sim_donation(switchtorho==1)=rho(switchtorho==1);
sim_donation(stay_default==1)=def(stay_default==1);
sim_donation(optout==1)=0;
add_zeros=zeros(addz, 1);
sim_donation=vertcat(sim_donation, add_zeros);


means=mean(sim_donation);
total=sum(sim_donation);


end